!MODULE PROXELMODULE
!The routine starts with initializing of some constants (e.g. CO2 in the atmosphere [ppm]), followed by the calls to the parameter-read-routines. The outputfile is opened (results.dat) and the Timestep-Loop is established. This loop runs from 1 to 20000 and a meteofileline is read in every loopstep. This means that the maximum number of lines in the meteofile is 20000 (but you can easily change this).
! In the loop there are the calls to GASFLUX and SOIL, and a line of output is written.


!crop***********************************************************************
      subroutine heats(jd,hr,tair,tb,sowdate,dtt)
!
!     calculates thermal time
!--------------------------------------------------------------------------

      USE Constants
      IMPLICIT NONE


      INTEGER :: jd,sowdate
      DOUBLE PRECISION :: hr,tair,tb,dtt,tav,tairsum

      save tairsum

      dtt=0.
      tav=0.

      if (int(hr+0.5).eq.0) then
         tairsum=tair
      else
         tairsum=tairsum+tair
      endif

      if (int(hr+0.5).eq.23) then
         tav=tairsum/24.
      endif

      if (jd.ge.sowdate) then
         if (tav.ge.tb) then
            dtt=tav-tb
         endif
      endif

      return
      end



!csoil--------------------------------------------------------------------
      subroutine proxel(styp,intercmax,bhoe, &
          intrv,xtair,wind,humidity,xsst,radsoil, &
          newLAI,gapf,atmpress,rain,trans,sresp,vint,xpsi1, &
          xpsi2,xpsi3,thet1,thet2,thet3,soiltmp1,soiltmp2,soiltmp3)

!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none

      double precision eta

      parameter (eta=0.6)

        double precision thet1,thet2,thet3,soiltmp1,soiltmp2,soiltmp3
      double precision bev2,intrv,klight,Rn_Gsoil,pyran
      integer k,styp
      double precision aerowid,bz0,bdd
      double precision thetaPsi,psiTheta,xtair,wind,humid,xsst,radsoil,  &
          newLAI,trans,sresp,humidity
      double precision atmpress,rain0,bra,bhoe,intercmax,gapf,  &
          thrf,vbes,vint,pEvSoil,rain,xpsi,xpsi1,xpsi2,xpsi3

      print*, 'in function proxel'


      vbes=0.
      vint=0.

!c     Globalstrahlung in W/m2
      pyran=radsoil*1000.

!c     Initialisieren
      bev2=-5.

!c     Falls Speicher voll
      if (interc.ge.intercmax) then
!c        Faellt alles durch
         thrf=rain
         interc=intercmax
!c     Falls Speicher nicht voll
      else
!c        Regen wird i.Abh. von gap-fraction interzeptiert
         interc=interc+(1.-gapf)*rain
!c        Speicher uebervoll?
         if (interc.ge.intercmax) then
            thrf=(interc-intercmax)+gapf*rain
            interc=intercmax
         else
            thrf=gapf*rain
         endif
      endif

!c     Bestandesniederschlag (rain0)
      rain0=thrf

      humid=humidity*0.01

!c     Verdunstung durch den Bestand
      vbes=trans

!c     Interzeptionsverdunstung
!c     aerodynamischer Widerstand 2m Hoehe
!c     Zo und Verdraengungshoehe berechnen
      if (bhoe.le.0.01) bhoe=0.01
      if (wind.le.0.01) wind=0.01
      if (bhoe.le.1.) then
         bz0=0.13*bhoe
      elseif (bhoe.gt.1) then
         bz0=(0.13-0.004*bhoe)*bhoe
      endif
      if (bz0.lt.0.08*bhoe) bz0=0.08*bhoe
      bdd=0.63*bhoe

!c     aerodynamischer Widerstand
      bra=aerowid(wind,bhoe,bz0,bdd)

!cMR   Hier Rn=-40.9+0.83pyran statt pyran verwenden
!cMR   und weiter Rn_G=Rn-hflux(1) (alles Wm-2 ???)
      call penmon(xtair,humid,-40.9+0.83*pyran,vint,intrv,  &
          atmpress,0.d0,bra)

!c     je nasser der Bestand, desto mehr naehert sich Interzept.verd.
!c     der maximal moeglichen (durch penmon berechneten)
!c     (die Konstante bev2 ist +/- willkuerlich & provisorisch)
!cMR?? woher kommt diese Gleichung, sie fuehrt dazu, dass interc
!cMR   nie mehr null wird ???
!cMR   soll wahrscheinlich Verdunstung aus unteren Schichten modellieren, aber
!cMR   so ganz geht das nicht, denn interc kann auch klein sein, wenn nur
!cMR   oberste Schicht gefuellt
      if (intercmax.gt.0.) then
         vint=vint*(1-exp(bev2*interc/intercmax))
      else
         vint=0.
      endif

!c     neuer Fuellzustand des Interzeptionsspeichers
      if (vint.le.interc) then
!c     interc wird um vint geleert
         interc=interc-vint
         if (interc.gt.intercmax) interc=intercmax
      else
         vint=interc
         interc=0.
      endif

!c     Computation of potential evaporation from soil
!c     zur Zeit noch voellig bulk
      if (bhoe.le.0.01) bhoe=0.01
      if (wind.le.0.01) wind=0.01
      if (bhoe.le.1.) then
         bz0=0.13*bhoe
      elseif (bhoe.gt.1) then
         bz0=(0.13-0.004*bhoe)*bhoe
      endif
      if (bz0.lt.0.08*bhoe) bz0=0.08*bhoe
      bdd=0.63*bhoe

      bra=aerowid(wind,bhoe,bz0,bdd)

!c     Empir. Bestimmung der verfuegbaren Energie zur Evap. Rn-G an der
!c     Bodenoberfl.
!c     extinkt.koef klight in ersten Term geschaetzt aus Puechabon-Daten
!c     LAI=2.9,
!c     median(i/i0)=0.19
!c     0.2 sei albedo von Boden
      klight=0.57
      Rn_GSoil=(1.-0.2)*pyran*exp(-klight*newLAI)-hflux(1)

!c     Aerodyn. Widerstand vom Boden
!c     einfache Annahme invers prop. zum windspeed, 600. ist empirisch
!c     so dass rd bei cpz um die 100s/m ist (vgl. Schelde 1996)

!cmr   Was kann man machen (wenn Bodenverdunstung zu hoch)?
!cmr     1) Faktor bhoe/15. rausnehmen (Faktor war 1 fuer Castelporziano), sollte
!cmr        vielleicht immer 1 sein. also 600./wind
!cmr     2) Wenn dass nichts hilft, kann man die 0. hochsetzen
!cmr        (Bodenwiderstand), am besten in Abhaengigkeit von Wassergehalt der
!cmr        obersten Schicht, z.B.
!cmr        rBoden=100./((ttheta-thetaR(i))/(thetaS(i)-thetaR(i))) aber vielleicht
!cmr        reicht es auch einfach auf 100. zu setzen (habe keine Erfahrung),
!cmr     3) Den Energiebilanzansatz von Eva einbauen.
!cgy      soilresistance=100./((theta(col,row,1)-thetaR(1))/
!cgy     +     (thetaS(1)-thetaR(1)))
      call penmon(xtair,humid,Rn_Gsoil,pEvSoil,intrv,atmpress,0.d0,   &
          600./wind)
!cgy     +     600.*bhoe/15./wind)

!c     Computation of unsaturated water flow mit rain-pEvSoil
!c     als oberer Randbedingung
      call unsatWaterFlow(rain0-pEvSoil,1.d0,styp)

!c     Teilverdunstungen summieren
      trans=vbes+vint+actEvap

!c     HEAT FLOW
      call compThermalProps

!c     call implicitHeatFlow(eta,3600.d0,xsst,20.d0)    ! 20.=deepTemp
      call implicitHeatFlow(eta,3600.d0,xsst)

!c     Root Uptake
      call qRootUp(vbes)

      do k=1,nOflayers
         theta(k)=theta(k)-RootUp(k)*1.0/   &
             (1000.0*thick(k))
      enddo

      do k=1,nofLayers
         if (theta(k).lt.(thetaR(k)+   &
              10.0**(-20.0*vgM(k)))) then
            write(*,*) 'Layer',k,' : Plants extract too much water'
            write(*,*) 'Theta dropped to',theta(k)
            write(*,*) 'PSI is now set to -10E5 ==> Mass Bal. not OK'
            psi(k)=-10000.
            theta(k)=ThetaPsi(psi(k),k)
         else
            psi(k)=psiTheta(theta(k),k)
         endif
      enddo

!c     Root water potential: PsiR=PsiSoil-q*(Rsoil+Rrootsurf)
!c     wobei Rrootsurf in anlehnung an Root Contact model berechnet
!c     rootSurfRes=thetaS/theta*(1/Rootdens)/permeability/thick
!c     cm-2                cm2cm-1root cmWS-1 h-1 cm

      do k=1,noflayers
         psiRoot(k)=psi(k)-rootUp(k)*   &
!c             Soil Resistance
             (1.0/(soilcond(k)*thick(k)*1000.)+  &
!c             Root surface resistance, root perm: 0.01 cm3 cm-1root cm-1WS
             thetaS(k)/(0.01*theta(k)*rootdens(k)*  &
             thick(k)*1000.))
      enddo

      call effectivePsi

!c     Berechnung der Bodenrespiration
!c     resptemp=temp(col,row,2)
      call compSoilResp(newLAI)
!cgy      sresp=0.0
!cgy      do k=1,nOflayers
!cgy         sresp=sresp+resp(k)       ! mmol m-2 h-1
!cgy      enddo
      sresp=resp(1)       ! mmol m-2 h-1

      thet1=theta(1)
      thet2=theta(2)
      thet3=theta(3)

      soiltmp1=temp(1)
      soiltmp2=temp(2)
      soiltmp3=temp(3)

!c     h2o-coupling; xpsi goes to report
      xpsi1=psi(1)
      xpsi2=psi(2)
      xpsi3=psi(3)

      return
      end





!csoil--------------------------------------------------------------------
      function psiTheta(ttheta,i)

!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none

      integer i
      double precision psiTheta,ttheta

!c     vgAlpha is [1/cm] ==> psi is [cm]
      psiTheta=(((ttheta-thetaR(i))/(thetaS(i)-thetaR(i)))  &
          **(-1.0/vgM(i))-1.0)**(1.0/vgN(i))/(-vgAlpha(i))

      return
      end





!csoil--------------------------------------------------------------------
      function kPsi(ppsi,i)

!c------------------------------------------------------------------------

      USE soil1pix_dble



      implicit none

        integer i
      double precision kPsi,ppsi,alphaPsi

!c     kSat in [mm/h] ==> kPsi in [mm/h]

!c     -Psi, da in van Genuchten Psi als pressure head (also positiv) ist
      alphaPsi=vgalpha(i)*(-ppsi)
      if (ppsi.gt.-1000.) then
         kPsi=ksat(i)*(1.0-alphaPsi**(vgN(i)-1.0)*   &
             (1.0+alphaPsi**(vgN(i)))**(-vgM(i)))**2.0/   &
             (1+(alphaPsi)**(vgN(i)))**(0.5*vgM(i))
      else
         kPsi=ksat(i)*(alphaPsi**(-0.5*vgM(i)*vgN(i))*   &
             (1-(1-alphaPsi**(-vgN(i)))**vgM(i))**2.0)
      endif

      return
      end





!csoil--------------------------------------------------------------------
      function thetaPsi(ppsi,i)

!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none


      integer i
      double precision thetaPsi,ppsi

!c     vgAlpha is [1/cm] ==> psi is [cm]
      thetaPsi=thetaR(i)+(thetaS(i)-thetaR(i))/  &
          ((1.0+(-vgAlpha(i)*ppsi)**vgN(i))**vgM(i))

      return
      end



!csoil--------------------------------------------------------------------
      function MaxPsiForDT(minDT)

!c     Berechnet das maximale Psi, welches bei einem gegebenem Boden
!c     und gegebenen minimalem Zeitschritt erreicht werden darf, um
!c     keine numerischen Instab. zu erzeugen
!c     Strategie ist also: modellierer gibt kleinsten Zeitschritt an, den
!c     er noch akzeptiert, und die Funktion legt dann den Wassergehalt fest
!c     bis zu dem jeder Layer aufgefuellt wird (Bucket-Modell-Konzept,
!c     ohne unbedingt Feldkapazitaet nehmen zu muessen
!c------------------------------------------------------------------------
      USE soil1pix_dble

    implicit none

        integer i
      double precision MaxPsiForDT,minDT,dt,maxPsiLayer,dz(nsl),   &
          ttheta,expo,maxPsi,dk_dt,dPsi_dT,kTheta

!c     Fuer jede Bodenschicht wird ausgerechnet ab welchem Matrixpotential
!c     der nach Moldrup et al. 1989 sich ergebende Zeitschritt kleiner
!c     wird als der von Benutzer vorgegebene minimale minDt
!c     dazu wird beim PWP angefangen, und dann PSI exponentiell (*0.8)
!c     erhoeht (PSI ist negetiv!)

      maxPsi=-(10.0**4.5)
      dt=minDt+1.0                                ! Stunden
      do i=1,noflayers-1
         dz(i)=(thick(i)+thick(i+1))/0.02         ! cm
         maxPsiLayer=-(10.0**4.5)     ! tiefer als PWP soll es aber nie fallen
 100     if ((dt.gt.minDT).and.(maxPsiLayer.lt.-1.0)) then
!c           vgAlpha is [1/cm] ==> maxpsiLayer is [cm]
            ttheta=thetaR(i)+(thetaS(i)-thetaR(i))/   &
                ((1.0+(-vgAlpha(i)*maxpsiLayer)**vgN(i))**vgM(i))
            expo=exp(dK_dT(ttheta,i)*dz(i)/   &
                (kTheta(ttheta,i)*dPsi_dT(ttheta,i)))
            dt=(expo-1.0)/(expo+1.0)*dz(i)/dK_dT(ttheta,i)
            maxPsiLayer=maxPsiLayer*0.8
            goto 100
         endif
         maxPsi=max(maxPsi,maxPsiLayer/0.8)
      enddo
!c     aus numerischen Gruenden soll Psi aber nie ueber -5.0cm
      MaxPsiForDt=min(maxPsi,-5.0)

      return
      end


!csoil--------------------------------------------------------------------
      function dPsi_dT(ttheta,i)

!c------------------------------------------------------------------------

      USE soil1pix_dble


      implicit none


      integer i
      double precision dPsi_dT,ttheta,dPsidT,effSat

      effSat=(ttheta-thetaR(i))/(thetaS(i)-thetaR(i))

      dPsidT=(effsat**(-1.0/vgM(i))-1.0)**(1.0/vgN(i)-1.0)*   &
          effSat**(-1.0/vgM(i)-1.0)/(-vgM(i)*vgN(i))
!c     -dPsi denn ohne "-" ist es dh/dT (van Genuchten pressure head)
      dPsi_dT=-dPsidT/(vgAlpha(i)*(thetaS(i)-thetaR(i)))

      return
      end





!csoil--------------------------------------------------------------------
      function kTheta(ttheta,i)

!c------------------------------------------------------------------------
      USE soil1pix_dble

    implicit none

        integer i
      double precision ttheta
      double precision kTheta,effSat

      effSat=(ttheta-thetaR(i))/(thetaS(i)-thetaR(i))
      kTheta=ksat(i)*sqrt(effSat)*(1.0-(1.0-     &
          effsat**(1.0/vgM(i)))**vgM(i))**2.0
      if (ktheta.eq.0.) then
         kTheta=1.e-15
         write(*,*) 'kTheta ',kTheta
      endif

      return
      end





!csoil--------------------------------------------------------------------
      function dk_dT(ttheta,i)

!c------------------------------------------------------------------------

      USE soil1pix_dble


      implicit none


      double precision dk_dT,ttheta
      double precision dKdT,effSat,inner,outer
      integer i

!c     Einheit wie kSat=mm/h
      effSat=(ttheta-thetaR(i))/(thetaS(i)-thetaR(i))

      inner=1.-effsat**(1./vgM(i))
      outer=(1.-inner**(vgM(i)))

      dKdT=0.5/sqrt(effSat)*outer**2+sqrt(effSat)*2.0*outer*   &
          inner**(vgM(i)-1.0)*effSat**(1.0/vgM(i)-1.0)

      dK_dT=dKDt/(thetaS(i)-thetaR(i))*kSat(i)
      if (dk_dt.eq.0.) write(*,*) 'dk_dt=',dkdt

      return
      end





!csoil--------------------------------------------------------------------
      subroutine penmon(bt,bfeuch,bstrge,evap,intrv,atmpress,rcxx,   &
          braxx)

!c     PENMAN-MONTEITH Verdunstung
!c------------------------------------------------------------------------
      implicit none

      double precision bt,bfeuch,bstrge,evap,intrv,atmpress
      double precision lhv,spw,psyko,pardif,dilu,rcxx,braxx,slope
      double precision tkel,wasmb,pardru

!c     SPW   = spez. Waerme Luft
!c     LHV   = latente Verdunstungswaerme Wasser
!c     PSYKO = Psychrometerkonstante
!c     PARDRU= aktuelles Saettigungsdefizit der Luft (hPa)
!c     DILU  = Dichte der Luft
!c     BRAXX = aerodyn. Widerstand
!c     RCXX  = Stomatawiderstand
!c     SLOPE = Aenderung des Saettigungsdampfdrucks

!c     Berechnung von
!c     spez. Waerme Wasser                 [J/kg*!]     (SPW)
!c     Verdampfungswaerme Wasser           [10e6*J/kg]  (LHV)
!c     Steigung der Saettigungskurve       [hPa/!]      (SLOPE)
!c     Psychrometerkonstante               [hPa/!]      (PSYKO)
!c     Wasserdampfsaettigungsdefizit,vpd   [hPa]        (PARDIF)
!c     akt. Wasserdampfpartialdruck        [hPa]        (PARDRU)
!c     Dichte der feuchten Luft                         (DILU)

!c     Temperatur in Kelvin
      tkel=bt+273.16

!c     spez. Waerme Wasser
      spw=16.92*exp(-0.6843*bt)+4171.54

!c     Anpassung an Luft [J/(kg*K)]
      spw=spw/4.197

!c     Verdampfungswaerme [J/kg]
      lhv=5.147*exp(-0.0004643*bt)-2.6466
      lhv=lhv*1000000.

!c     Steigung des Saettigungsdampfdruckes
      wasmb=6.1078*exp((17.08085*bt)/(234.175+bt))

!c     slope of vapour pressure curve [hPa/K]
      slope=0.6341*exp(0.04683*bt)-0.17268

!c     VPD
      pardif=wasmb*(1.-bfeuch)
!c     Wasserdampfpartialdruck
      pardru=wasmb*bfeuch

!c     Dichte der feuchten Luft [kg/m^3]
      dilu=(atmpress/(287.04*tkel))*(1-((0.378*pardru)/atmpress))
      dilu=dilu*100.

!c     Psychrometerkonstante
!c     from Monteith/Unsworth [hPa/K]
      psyko=spw*atmpress/(lhv*.622)

!c     PENMAN-MONTEITH formal berechnen
      if (braxx.eq.0) then
         write(*,*) 'PENMON: braxx.eq.0.'
         stop 'flush the toilet'
      endif

!c     Penman-Monteith [kg/(m^2*s)]
      evap=(slope*bstrge+spw*dilu*pardif*(1./braxx))/  &
          (slope+psyko*(1.+(rcxx/braxx)))              ! Jm-2s-1
      evap=evap/lhv

!c     Umrechnung in mm/Zeitintervall
      evap=evap*intrv

      return
      end





!csoil--------------------------------------------------------------------
      function aerowid(uu,hh,z0,dd)
!c     aerodynamic resistance [s/m]
!c     from Mauser's original PointModel
!c------------------------------------------------------------------------
      implicit none

      double precision vonka2,ra1
      double precision mhelp
      double precision aerowid,uu,z0,dd,hh

      vonka2=0.16

!c     aerodynamischer Widerstand
      mhelp=((hh+2.)-dd)/z0
      if (mhelp.le.0.0) then
         write(*,*) 'Aerowid: mhelp < 0.'
         read(*,*)
         aerowid=1.e6
         return
      endif
      ra1=(real(log(mhelp))**2)/(vonka2*uu)
      aerowid=ra1

      return
      end






!csoil--------------------------------------------------------------------
      subroutine unsatWaterFlow(q0pot,global_dt,styp)

!c     Computes unsaturated water flow by moving mean slope model
!c     of MOLDRUP et al. 1989
!c------------------------------------------------------------------------

      USE soil1pix_dble


      implicit none


      double precision psiMin

      parameter (psiMin=-1.0d10)

      integer i,cracks,styp
      double precision alphaL(nsl),kL(nsl)
      double precision alphaN(nsl),kN(nsl),dz(nsl)
      double precision qInLayAkt,dpsi,psi0Max,  &
          qMax,q0akt,q0pot,q0potplus,psimax,   &
          alphadeep,kLdeep
      double precision dt,elapsTime,global_dt
      double precision thetaPsi,kPsi,dk_dt,psiTheta

      cracks=0
!c      pond=0.0
      runoff=0.
      actevap=0.0
      dpsi=0.2
!c      dTheta=0.005
      psimax=-5.
      elapsTime=0.0
 100  if (elapstime.lt.global_dt) then
         do i=1,nOfLayers
!c           eq 10, but for velocity reasons dpsi set to 0.5 X
!c           psi, dpsi in [cm] ==> alpha [1/cm]
!c           k=kPsi(psi(i),i)
            if (psi(i).le.psiMin) then
               psi(i)=psiMin
               theta(i)=thetaPsi(psiMin,i)
               write(*,*)'Soil Psi below -10e10: adjusted'
            endif
            if ((psi(i)-dpsi).gt.-1000.) then
               if (kPsi(psi(i)+dpsi,i).le.0.0)  &
                   write(*,*) 'alphaL(i)-log-prob.'
               alphaL(i)=log(kPsi(psi(i)+dpsi,i)/   &
                   kPsi(psi(i)-dpsi,i))/(2.0*dpsi)
            else
               alphaL(i)=-(3.*vgn(i)-1.0)/psi(i)
            endif
!c           kPsi in [mm/h]
            kL(i)=kPsi(psi(i),i)/  &
                exp(alphaL(i)*psi(i))
         enddo

!c        Computation of qH2O unsaturated for layer 1 to n-1
         do i=1,nOfLayers-1
            alphaN(i)=(alphaL(i)+alphaL(i+1))/2.0
            kN(i)=(kL(i)+kL(i+1))/2.0

!c           Wenn Matrix-flux potential concept benutzt werden soll,
!c           dann korrektur einfuegen
!c           vgl. Shaykewich et al. 1977, eq. 13
!c           if (psi(i).ne.psi(i+1)) then
!c              k=kn(i)/(0.5*alphaN(i)*(psi(i)-psi(i+1))*
!c     +           (1+2*exp(alphaN(i)*psi(i+1))/
!c     +           (exp(alphaN(i)*psi(i))-
!c     +           exp(alphaN(i)*psi(i+1)))))
!c              write(*,*) kn(i),k
!c           endif
            dz(i)=(thick(i)+thick(i+1))/0.02            ! cm
!c           kN [mm/h] ==> qH2O is in [mm/h]

            qH2O(i)=kN(i)*(exp(alphaN(i)*psi(i))-   &
                exp(alphaN(i)*psi(i+1)))/   &
                (exp(alphaN(i)*dz(i))-1)+   &
                kN(i)*exp(alphaN(i)*psi(i))
         enddo


!c        Computation of qH2O unsaturated for deepest layer
         if (lowbouZero.gt.0) then
!c           Zero head boundary condition:
!c           Es ist eine Schicht drunter mit gleicher hydr.
!c           Charakteristik aber psimax-Potential
            i=noflayers
            dz(i)=thick(i)/0.02     ! Groundwater directly on Unterer Grenze
            alphaDeep=log(kPsi(psimax+dpsi,i)/  &
                kPsi(psimax-dpsi,i))/(2.0*dpsi)
!c           kPsi in [mm/h]
            kLdeep=kPsi(psimax,i)/exp(alphaDeep*psimax)
            alphaN(i)=(alphaL(i)+alphaDeep)/2.0
            kN(i)=(kL(i)+kLdeep)/2.0

            qH2O(i)=kN(i)*(exp(alphaN(i)*psi(i))-   &
                exp(alphaN(i)*psimax))/   &
                (exp(alphaN(i)*dz(i))-1)+  &
                kN(i)*exp(alphaN(i)*psi(i))

!c           qh2o(i)=-1.0
         else
!c           Unit gradient boundary condition
!c           lower layer: kN=kL, da nach unten unit gradient free
!c           drainage Boundary condition
            qH2O(nOfLayers)=kL(nOfLayers)*exp(alphaL(noFLayers)*   &
                psi(nOfLayers))
         endif

!c        groesser darf dt nicht werden
         dt=global_dt-elapsTime

         do i=1,noflayers-1
!c           for each layer (vektor)
!c           dz is [cm], kTheta [mm/h] ==> dt is in [h]
!c           (10 mm/cm *2.0=factor 20.0)

!c           dt=min(dt,dz(i)*(exp(alphaN(i)*dz(i))-1)/(exp(alphaN(i)*
!c     +        dz(i))+1)*20.0*dTheta/(kTheta(theta(col,row,i)+dTheta,i)-
!c     +        kTheta(theta(col,row,i)-dTheta,i)))
!c           fuer korrekte Steuerung bei Infiltration vielleicht noch
!c           ein alpha fuer infiltr. in oberste bodenschicht berechnen!
            dt=min(dt,dz(i)*(exp(alphaN(i)*dz(i))-1)/  &
                (exp(alphaN(i)*dz(i))+1)/  &
                dK_dT(theta(i),i))
!cdebug
            if (dt.le.0.) then
               write(*,*) 'dt-error',i,elapsTime,dt,   &
                   dK_dT(theta(i),i),theta(i),   &
                   dz(i),alphan(i)
!cwatz               read(*,*)
            endif
         enddo

!c        if ((q0pot.gt.0.0).and.(thetaMax(1).gt.thetaS(1)*0.9)) then
!c           write(*,*) dt
!c           Wenn Infiltration, dann nur so, dass Boden nicht staerker
!c           als bis thetaS auf gefuellt wird
!c           dt=min(dt,max(1000.0*thick(1)/
!c     +        (ksat(1)-qh2o(1))*(thetaS(1)-theta(col,row,1)),0.001))
!c           write(*,*) dt
!c           write(*,*) '----------------'
!c           read(*,*)
!c        endif

!c        aufpassen dass theta nicht unter null sinkt??

!c        Computation of actual infiltration or Evaporation
         if ((q0pot.lt.0.0).and.pond.lt.0.1) then
!c           q0pot<0 : Evaporation from soil surface
!c           Computation of evaporation capacity of soil (nach N. Jarvis 1995)
            qSatFlow(1)=0.0
            psi0Max=((-thick(1)*bb(1)/2.0)-psi(1)*(1.0-bb(1)))/   &
                (1.0+bb(1))
            qMax=kSat(1)*(-0.5*vgAlpha(1)*(psi0Max+  &
                psi(1)))**bb(1)*  &
                ((psi(1)-psi0Max)/(5.0*thick(1))-1.0)
                                  ! 5.0* da faktor 10 fuer m in cm
!c            write(*,*) psi0Max,qMax,q0pot
!c            read(*,*)
            q0akt=max(-qmax,q0pot)         ! Bsp. max(-0.01,-0.5)=-0.01

            theta(1)=theta(1)+(q0akt-qH2O(1))*dt/  &
                (1000.0*thick(1))

            actEvap=actEvap-q0akt*dt
            infill=0.0
         else
!c           q0pot<=>0 and ponding water : Infiltration due to throughfall
            q0potplus=q0pot+pond/dt
            q0akt=min(ksat(1),q0potplus)

            infill=q0akt
            actEvap=0.0
            if (cracks.eq.0) then
               pond=pond+(q0Pot-q0akt)*dt
               runoff=runoff+(q0Potplus-q0akt)*dt
            endif
!c           wenn Runoff nicht rein soll, dann q0akt

            theta(1)=theta(1)+(q0akt-qH2O(1))*   &
                dt/(1000.0*thick(1))

            if (theta(1).gt.thetaMax(1,styp)) then
               qSatFlow(1)=(theta(1)-thetaMax(1,styp))*   &
                   (1000.0*thick(1))/dt
               theta(1)=thetaMax(1,styp)
            else
               qSatFlow(1)=0.0
            endif
!c           simple Simulation of Macropore flow
            if (cracks.eq.1) then
               qsatflow(1)=qsatflow(1)+(q0Potplus-q0akt)
!c              i.e. all what is above Ksat infiltrates into next layer!
            endif
         endif
         psi(1)=psiTheta(theta(1),1)

!c        Updating theta of deeper layers by explicit integration
!c        and by allowing for fast unidirectional drainage
         do i=2,noflayers
!c           k=theta(col,row,i)
            qInLayakt=min(ksat(i),qSatFlow(i-1))         ! if runoff
!c           qinLayakt=qSatFlow(i-1)
            runoff=runoff+(qSatFlow(i-1)-qInLayakt)  &
                *dt

            infill=infill+qinLayakt
!c           if (qSatFlow(i).gt.0.0)
!c     +        write(*,*) 'qs,qinakt, run',qSatFlow(i),qinLayakt,runoff(col,row)
            theta(i)=theta(i)+(qH2O(i-1)-qH2O(i)+  &
                qInLayAkt)*dt/(1000.0*thick(i))
            if (theta(i).gt.thetaMax(i,styp)) then
               qSatFlow(i)=(theta(i)-thetaMax(i,styp))*  &
                   (1000.0*thick(i))/dt
               theta(i)=thetaMax(i,styp)
            else
               qSatFlow(i)=0.0
            endif
            if (cracks.eq.1) then
               qsatflow(i)=qsatflow(i)+(qSatFlow(i-1)-qInLayakt)
               runoff=0.0
            endif
            psi(i)=psiTheta(theta(i),i)
         enddo

         if (cracks.eq.1) then
            pond=max(pond+(q0pot-infill)*dt,0.0)
         else
            pond=pond-runoff
         endif

!cgy???
         drainage=(qH2O(noflayers)+qSatFlow(noflayers))*dt

         elapstime=elapstime+dt
!c        write(*,*) 'elaps',elapstime,'dt',dt
!c        read(*,*)
         goto 100
      endif

      return
      end






!csoil--------------------------------------------------------------------
      subroutine compThermalProps

!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none


      integer i
!c     double precision  a,b,c,d,e,volFracSolid
!c     double precision  clay


!c     clay=0.12

      do i=1,nOfLayers
!c        volFracSolid=bd(i)/2.65

!c        Heat Capacity
         spHeatCap(i)=2.4E+6*(1.-thetaS(i)/0.96)+4.18E+6*   &
             theta(i)

!!c        Thermal Conductivity (vereinfacht nach Campbell 1985, S.33ff
!!c        a=0.65-0.78*BD(i)+0.60*bd(i)*BD(i)
!c        b=2.8*theta(col,row,i)*volFracSolid
!c        c=1.0+2.6/sqrt(clay)
!c        d=0.03+0.7*volFracSolid*volFracSolid
!c        e=4.0

         lambTemp(i)=lambA(i)+lambB(i)*theta(i)*  &
             theta(i)-(lambA(i)-lambD(i))*  &
             exp(-(lambC(i)*theta(i))**lambE(i))

!c        write(*,'(i3, 4F12.2)') i,bd(i),theta(col,row,i),spHeatCap(i),
!c    +      lambTemp(i)
      enddo

      return
      end



!csoil--------------------------------------------------------------------
!c      subroutine implicitHeatFlow(col,row,eta,dt,uppBoun,deepTemp)
      subroutine implicitHeatFlow(eta,dt,uppBoun)

!c------------------------------------------------------------------------
      USE soil1pix_dble

      implicit none

      integer i
      double precision eta,einsMinEta,dt,avgLambda,uppBoun,Tsurf,k0
      double precision e(nsl),f(nsl),d(nsl),r(nsl)
      double precision T(nsl),Tnew(nsl),k(nsl),heatCap(nsl)


      einsMinEta=1.-eta

      do i=1,nOflayers
         T(i)=temp(i)
      enddo

      avgLambda=(lambTemp(1)*thick(1)+lambTemp(2)*0.5*thick(2))/  &
          (thick(1)+0.5*thick(2))
      k(1)=avgLambda/(zz(2)-zz(1))
      heatCap(1)=spHeatCap(1)*thick(1)/dt

!c     Neumann condition (flux into first layer given)
!c     e(1)=-eta*k(1)-heatCap(1)
!c     f(1)=eta*k(1)
!c     r(1)=T(1)*(-heatCap(1)+einsminEta*k(1))-k(1)*einsMinEta*T(2)-
!c    +     uppBoun

!c     Dirichlet condition (Temperature first layer given)
!c     e(1)=1
!c     f(1)=0
!c     r(1)=uppBoun              ! uppBoun hier als Temp interpretiert

!c     New approach: Tsurf bekannt
      tSurf=uppBoun
      k0=lambtemp(1)/(thick(1)*0.5)
      e(1)=-eta*(k0+k(1))-heatCap(1)
      f(1)=eta*k(1)
      r(1)=T(1)*(einsmineta*k0+einsmineta*k(1)-heatCap(1))-k0*Tsurf-   &
          k(1)*T(2)*einsmineta

      do i=2,nOfLayers-1
         heatCap(i)=spHeatCap(i)*thick(i)/dt
         avgLambda=(lambTemp(i)*thick(i)+lambTemp(i+1)*thick(i+1))/   &
             (thick(i)+thick(i+1))
         k(i)=avgLambda/(zz(i+1)-zz(i))

         f(i)= eta*k(i)
         d(i)=f(i-1)
!c        Wenn Dirichlet condition dann
         d(i)=eta*k(i-1)
         e(i)=-eta*(k(i)+k(i-1))-heatCap(i)
         r(i)=-einsMinEta*k(i-1)*T(i-1)-(heatCap(i)-einsMinEta*(k(i)+   &
             k(i-1)))*T(i)-einsMinEta*k(i)*T(i+1)
      enddo

      heatCap(nofLayers)=spHeatCap(nofLayers)*thick(nofLayers)/dt
      k(nOfLayers)=lambTemp(nofLayers)/thick(nofLayers)

!c     Dirichlet condition
!c     d(nOfLayers)=0
!c     e(nOfLayers)=1
!c     r(nOfLayers)=deeptemp

!c     Neumann cond. with no flow at lower boundary
      d(nOfLayers)=eta*k(nOfLayers-1)
      e(nOfLayers)=-eta*k(nOfLayers-1)-heatCap(nOfLayers)
      r(nOfLayers)=0.0-heatCap(nOfLayers)*T(nOfLayers)+einsMinEta*   &
          k(nOfLayers-1)*(T(nOfLayers)-T(nOfLayers-1))


!c     Thomas algorithmus
      do i=1,nOfLayers-1
         f(i)=f(i)/e(i)
         r(i)=r(i)/e(i)
         e(i+1)=e(i+1)-d(i+1)*f(i)
         r(i+1)=r(i+1)-d(i+1)*r(i)
      enddo

      Tnew(nOfLayers)=r(nOfLayers)/e(nOfLayers)
      do i=nOfLayers-1,1,-1
         Tnew(i)=r(i)-f(i)*Tnew(i+1)
      enddo
      do i=1,nOfLayers-2
!c        hFlux(i) wird jetzt als Fluss IN den Layer i interpretiert!
         hFlux(i+1)=-k(i)*(eta*Tnew(i+1)+(1-eta)*T(i+1)-eta*Tnew(i)-   &
             (1-eta)*T(i))
!c        write(*,'(i3,5f10.4),/') i,heatCap(i),k(i),T(i),Tnew(i),hFlux(i)
      enddo
      hFlux(1)=-k0*(eta*Tnew(1)+(1-eta)*T(1)-eta*Tsurf-(1-eta)*Tsurf)
      i=nOfLayers
      hFlux(i)=0.0
!c     write(*,'(i3,5f10.4),/') i,heatCap(i),k(i),T(i),Tnew(i),hFlux(i)

!c     write(*,*) '--------------------------------------'
      do i=1,nOflayers
         Temp(i)=Tnew(i)
      enddo

      return
      end


!csoil--------------------------------------------------------------------
      subroutine qRootUp(transp)

!c     Berechnung der root-uptake nach Layer, wobei Transpiration
!c     gegeben ist.
!c     Ansatz: qRoot(i)=k(Psi(i))*(PsiSoil(i)-PsiPlant) und
!c     summe(qRoot)=Transp da PsiPlant sehr gross, wird PsiSoil-PsiPlant
!c     zu konst. deltaPsi zusammengef.,wonach die Gleichungsys. ohne
!c     negative Fluesse oder dgl. loesbar istund sich eine mit der
!c     soilCond. gewichtete Wasseraufnahme ergibt.
!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none

      integer i,flag,nOfNewEmpty
      logical empty(nsl),oldempty(nsl)
      double precision qRootMax(nsl),qRootDemand(nsl),qRootakt(nsl)
      double precision transp,tRest,sumqrootdemand
      double precision soilRootCond
      double precision conduct(nsl),totalRes(nsl),sumconduct

!c     Volumen-Conductivity in Flaechenbezogene Conducti. umrechnen
!c     (cm*h)^-1 wird zu mm/(cm*h), d.h. nach Multiplik. mit
!c     Potentialdifferenzergibt sich statt cm3/cm3*h eben mm/h
!c     (Potential und cm in conduct kuerzen sich)

      do i=1,noflayers
         soilcond(i)=max(soilrootcond(i),1.0e-25)

         totalRes(i)=   &
!c             Soil Resistance mm cm-1Ws h-1
             1.0/(soilcond(i)*thick(i)*1000.)+   &
!c             Root surface resistance, root perm: 0.01 cm3 cm-1root cm-1WS
             thetaS(i)/   &
             (0.0001*theta(i)*rootdens(i)*thick(i)*1000.)

         conduct(i)=1./totalRes(i)
!c         write(*,*) i,conduct(i),transp,totalRes(i)

!c        alternativer Ansatz, Wasseraufnahme zu begrenzen: vergleichbarer
!c        Ansatz wie bei der Evaporation (vgl. N. Jarvis 1995)

!c        psi0Max(i)=psi(i)*(bb(i)-1.0)/(1.0+bb(i))
!c         r2(i)=1.0/sqrt(rootDens(i)*3.14159)
!c         r1(i)=min(0.1+0.0*r2(i),0.9*r2(i))
!c        cm Wurzelradius (+0.0*r2,um Vektor zu erhalten)

!c        cmThick*cmWS*cmRoots/cm3soil*mmH2O/h=mmH2O/h,
!c        daher thick(i)*100.0 (m --> cm)
!c        Vielleicht kann man hier mal Matrix-Flux Potential konzept
!c        einsetzen um bessere Schaetzungen zu erhalten, statt einfach
!c        Mittelwert der Konductivities zu Nehmen
!c        qRootMax(i)=2.0*314.15*thick(i)*rootDens(i)*kSat(i)*
!c    +      (-0.5*vgAlpha(i)*(psi0Max(i)+psi(i)))**bb(i)*
!c    +      ((psi(i)-psi0Max(i))/log(r2(i)/r1(i)))

!c        qRootEnd(i)=100.0*thick(i)*rootDens(i)*kSat(i)*
!c    +      (-0.5*vgAlpha(i)*(-1e+5+psi(i)))**bb(i)*
!c           ((psi(i)-(-1e+5))/log(r2(i)/r1(i)))

         qrootmax(i)=(theta(i)-thetaR(i)-0.005)*thick(i)*1000.
      enddo

      sumConduct=0.0
      do i=1,noflayers
         sumConduct=sumConduct+conduct(i)
      enddo

      if (sumConduct.gt.0) then
         do i=1,noflayers
            qRootDemand(i)=transp*(conduct(i)/sumConduct)
         enddo
      else
         write(*,*) 'sumConduct =< 0',sumConduct,transp,   &
             theta(1),soilcond
         read(*,*)
      endif

!c     Algorithmus um verbleibende Wasseraufnahme auf andere Schichten zu
!c     verteilen, wenn eine (oder mehrere) Schichten nicht so viel Wasser
!c     bieten koennen wie obige Zeile berechnet
      do i=1,noflayers
         empty(i)=.false.          ! Vektor ueber Bodenschichten, ob leer
         oldempty(i)=.false.
         if (transp.lt.1.0e-4) then
            trest=0.0
            qRootakt(i)=min(qrootDemand(i),qrootmax(i))
         endif
      enddo

      tRest=transp     ! restliche Transpiration, die noch zu verteilen ist
 100  if (tRest.gt.0.00001) then
         do i=1,noflayers
            if (.not.empty(i)) qRootakt(i)=min(qrootDemand(i),   &
                qrootmax(i))
            empty(i)=(qrootmax(i).lt.qrootakt(i))
         enddo

!c        Wenn keine neue Bodenschicht leer geworden ist
         nOfNewEmpty=0
         do i=1,noflayers
            if (empty(i).neqv.oldempty(i)) then
               nOfNewEmpty=nOfNewEmpty+1
               if (empty(i)) tRest=tRest-qRootakt(i)
               write(*,*) qrootdemand(i),qrootakt(i)
            endif
            if (nOfNewEmpty.eq.0) trest=0
            if (empty(i)) qRootDemand(i)=0.
         enddo

         sumqrootdemand=0.0
         do i=1,noflayers
            sumqrootdemand=sumqrootdemand+qrootdemand(i)
         enddo

         if (sumqrootdemand.gt.0.) then
            do i=1,noflayers
               qRootDemand(i)=qRootDemand(i)*tRest/sumqRootDemand
            enddo
         endif

         flag=0
         do i=1,noflayers
            if (.not.empty(i)) flag=1
         enddo
         if (flag.eq.0) goto 200

         do i=1,noflayers
            oldempty(i)=empty(i)
         enddo

         goto 100
      endif
 200  continue

      do i=1,noflayers
         RootUp(i)=qrootakt(i)
!c        Root water potential
!c        psiRoot(i)=psi(i)-qrootakt(i)/conduct(i)
      enddo

      return
      end






!csoil--------------------------------------------------------------------
      function soilRootCond(i)

!c     Berechnung der 1/(Soil Resistance) pro Volumen nach Moldrup et al.
!c     1992 eq 9 (no steady state assumption)
!c------------------------------------------------------------------------
      USE soil1pix_dble

    implicit none

      double precision psiMin

      parameter (psiMin=-1.0d10)

        integer i
      double precision soilrootcond,thetaorg,ttt1,ttt2,b,r1,r2,rs
      double precision psiTheta
      double precision kTheta

      thetaorg=theta(i)

      theta(i)=thetaorg+0.0005
      ttt1=PsiTheta(theta(i),i)
      theta(i)=thetaorg-0.0005
      ttt2=PsiTheta(theta(i),i)
      b=abs(log(ttt1/ttt2))/log((thetaorg+0.0005)/  &
          (thetaorg-0.0005))
      theta(i)=thetaorg

      if ((rootDens(i).eq.0.0).or.(psi(i).le.psiMin)) then
         soilrootcond=0.0
      else
         r2=1./sqrt(rootDens(i)*3.1415)
         r1=0.1                                        ! cm Wurzelradius
         rs=1.5/(b*0.31415926*ktheta(theta(i),i)*rootDens(i))*   &
             (log(r2/r1)-0.5+(r1*r1)/(r2*r2-r1*r1)*   &
             log(r2/r1))

!c        soilrootcond=0.66666666667*(b*0.31415926*ktheta(theta(col,row,i),i)*
!c    +      rootDens(i))/
!c    +      (log(1/(sqrt(rootDens(i)*3.1415)*r1))-0.5)
         soilrootcond=1./rs
!c        write(*,*) soilrootcond
!c        soilrootcond is in 1/(cm*h), ktheta in mm/h daher faktor 0.1*PI
      endif

      return
      end






!csoil--------------------------------------------------------------------
      subroutine effectivePsi

!c     Calculation of weighed soil Psi: weights are qRootup(i)/(qRootup(i)+b)
!c     where bDilu is the Tardieu (1993) dilution Parameter
!c------------------------------------------------------------------------

      USE soil1pix_dble

      implicit none


      integer i
      double precision bDilu,totUp,totRoots,effPsi1,effPsi2,effPsiTard
      double precision weights(nsl),sumweights


      bDilu=0.001
      totUp=0.0
      totRoots=0.0
      do i=1,noflayers
         totUp=totUp+rootUp(i)
         totRoots=totRoots+rootdens(i)*thick(i)
         weights(i)=rootUp(i)/(bDilu+rootUp(i))
      enddo

      effPsi1=0.
      effPsi2=0.
      effPsiTard=0.0
      sumweights=0.0
      do i=1,noflayers
         effpsi1=effpsi1+psiroot(i)*rootup(i)
         effpsi2=effpsi2+psiroot(i)*rootdens(i)*thick(i)
         effPsiTard=effPsiTard+psi(i)*weights(i)
         sumweights=sumweights+weights(i)
      enddo

      if (totup.gt.0.00001) then
!c        effpsi=0.5*(effpsi1/totup+effpsi2/totRoots)
         effPsi=effPsiTard/sumweights
      else
         effpsi=-9999.
      endif

      return
      end





!csoil--------------------------------------------------------------------
      subroutine compSoilResp(LAI)

!c     computeSoilRespiration
!c------------------------------------------------------------------------
      USE soil1pix_dble

    implicit none

      double precision mol_c

      parameter (mol_c=0.012011)

        integer i
      double precision ThetaSoil,E0,T0,tRef,thetaHalf,theta0,tempfac,   &
          H2Ofac,respNorm,LAI

      do i=1,nOfLayers
         thetaSoil=theta(i)
         E0=respTP1(i)
         T0=respTP2(i)
         tRef=respTP3(i)
         thetaHalf=respFP1(i)
         theta0=respFP2(i)

!c        Computation of environmental conditions for soil resp.
         tempfac=exp(E0*(1/(tref-t0)-1/(temp(i)-t0)))
         H2OFac=(thetaSoil-theta0)/((thetaHalf-Theta0)+   &
             (thetaSoil-theta0))
         H2Ofac=max(H2Ofac,0.01)
         envfac(i)=tempFac*H2OFac

!c        Normierte Atmung respNorm bei optimaler Feuchte und Tref
!c        corg is g/m3, ki is 1/h, thick is m: so resp is g/m2h
!c         respNorm=thick(i)*(corg(i)*k1(i)*fractK1(i)+
!c     +        corg(i)*k2(i)*(1.0-fractK1(i)))
!c        write(*,'(5f14.5)') corg(i),k1(i),thick(i),envfac(i),resp(i)

!c        solange kein Boden-Corg Info dieser Wert beste Schaetzung
!c        (2.3 umol m-2 s-1 zu 0.36 g m-2 h-1)
!c        untere Schichten sollen nicht atmen
         if (i.eq.1) then
            respNorm=2.3*mol_c*3.6     ! 2.3 hoeher -> resp hoeher
         else
            respNorm=0.
         endif

         resp(i)=envfac(i)*respNorm
      enddo

!c     Umrechnung von g m-2 h-1 in mmol m-2 h-1: MCO2 in kg/mol oder g/mmol
!c     Vektor durch Skalar
      do i=1,nOfLayers
         resp(i)=resp(i)/mol_c
      enddo

!cmr   --------------------------------------------------------------
!c     Bezueglich der Bodenatmung schlage ich vor, das boden.atm wegzulassen,
!c     und dafuer zunaechst folgende Formel anzugeben:
!c     Rsoil=(0.6+1.29*LAI)*theta(1)/(theta(1)+0.04)*exp(58.2+285.4*
!c           theta(1)/0.2*(1.0/(18.+46)-1.0/(Tsoil(1)+46)),
!c     wobei theta(1) Wassergehalt in erster Schicht ist und Tsoil entsprechend.
!c     Rsoil ist dann in umol m-2 s-1.
      resp(1)=(0.6+1.29*LAI)*theta(1)/(theta(1)+0.04)*   &
          exp((58.2+285.4*theta(1)/0.2)*(1.0/(18.+46)-1.0/   &
          (temp(1)+46)))
      resp(1)=resp(1)*3.6   ! umol m-2 s-1 -> mmol m-2 h-1
!c     --------------------------------------------------------------

      return
      end




!cdoc---------------------------------------------------------------------------
!cdoc                    Variablenliste PROXEL
!cdoc
!cdoc  VarName           Inhalt, Einheiten
!cdoc  -------------------------------------------------------------------------
!cdoc  atag              Anfangstag
!cdoc  atmpress          Luftdruck, [hPa]
!cdoc  bdd               zero plane displacement = 0.63*bhoe(ki) [m]
!cdoc                    Rauhigkeitslaenge der Landschaft
!cdoc  beva              Evaporation nach Penman-Monteith [mm/intrv]
!cdoc  bfeuch            Luftfeuchte [%]
!cdoc  bhoe              Vegetationshoehe [m]
!cdoc  bodart            Bodenart je nach ./bod/file
!cdoc                    z.B. 1=Ton,2=Tonlehm,3=Lehm,4=Lehmsand,5=Sandlehm,
!cdoc                         6=Sand,7=Moor,8=Fels
!cdoc  cpsi(l)           Saugspannung [MPa]
!cdoc  bra               aerodyn. Widerstand [s/m]
!cdoc  braxx             aerodyn. Widerstand [s/m]
!cdoc  bz0               roughness length [m]
!cdoc  ca                CO2-Konzentration ueber Bestand [ppm]
!cdoc  dilu              Dichte der Luft [kg/m3]
!cdoc  finterc           Faktor fuer maximale Interzeption intercmax [mm]
!cdoc  gapf              gap fraction, nicht von Pflanzen bedeckte
!cdoc                       Grundflaeche [-],
!cdoc                       = 1-(Kronenflaeche/Grundflaeche)
!cdoc  gfac              Skalierungsfaktor fuer stomataeren Leitwert [-]
!cdoc  gfacpot           maximaler gfac wenn Bodenwasser
!cdoc                      uneingeschraenkt verfuegbar
!cdoc  gfk1-6            Faktor fuer gfac als Funktion von cpsi
!cdoc                      (Bedeutung je nach gfv)
!cdoc  gfv               Variante der Beziehung zwischen gfac und cpsi:
!cdoc                          0=konstant gfk1
!cdoc                          1=negativ-logistisch (fest),
!cdoc                          2=linear (fest),
!cdoc                          3=10^irgendwas (fest),
!cdoc                          4=negativ-logistisch
!cdoc                          5=linear
!cdoc                          6=PowerCurve
!cdoc                          7=2 Geraden (3 Punkte)
!cdoc                          sonst keine (dann gfac=gfacpot)
!cdoc  interc            Interception, [mm]
!cdoc  intercmax         maximale Interzeption, [mm]
!cdoc  intrv             Zeitschritt in Sekunden (3600 sec)
!cdoc  lhv               latente Verdunstungswaerme Wasser [J/kg]
!cdoc  mhelp             Hilfsvariable aerodynamischer Widerstand
!cdoc  pardif            aktuelles Saettigungsdefizit der Luft (=vpd)
!cdoc                      (hPa) (heisst z.T. faelschlicherweise pardru)
!cdoc  pardru            akt. Dampfdruck [hPa]
!cdoc  pEvSoil           potentielle Evaporation vom Boden (mm/h)
!cdoc  psyko             Psychrometerkonstante [hPa/K]
!cdoc  rain              Freilandniederschlag [mm/h]
!cdoc  rain0             Freilandniederschlag, [mm/intrv]
!cdoc  rcxx              Oberflaechenwiderstand [s/m]
!cdoc  rst               stomataerer Leitwert oder (am Schluss)
!cdoc                       Widerstand [s/m]
!cdoc  slope             Aenderung des Saettigungsdampfdrucks mit der
!cdoc                      Temp [hPa/K]
!cdoc  spw               spez. Waerme Luft [J/(kg*K)]
!cdoc  thrf              Regen ohne Kronenkontakt auf Boden [mm/intrv]
!cdoc  tkel              Lufttemperatur [K]
!cdoc  vbes              Verdunstung durch den Bestand, [mm/intrv]
!cdoc  vint              Interzeptionsverdunstung, [mm/intrv]
!cdoc  vonka2            (von Karman-Konstante)**2 = 0.16
!cdoc  wasmb             Saettigungsdampfdruck [hPa]
!cdoc---------------------------------------------------------------------------
!cdoc  theta             Bodenwassergehalt   m3/m3
!cdoc  psi               Bodenmatrixpotential Saugspannung cmWS
!cdoc  qH20              ungesaettigtes Wasserfluss mm/h oder MPa
!cdoc  rootUp            Wasseraufnahme durch Wurzeln mm/h
!cdoc  SoilCond          Bodenleitfaehigkeit
!cdoc  psiRoot           Saugspannung der Wurzeln
!cdoc  temp              Bodentemperatur
!cdoc  hFlux             Heatflux Bodenwaermestrom
!cdoc  envfac            Skalar fuer Bodenatmung
!cdoc  resp              Bodenatmung pro Bodenschicht
!cdoc  rsoil             Gesamtbodenatmung
!cdoc  infill            Infiltration durch Regen
!cdoc  runoff            Oberflaechenabfluss
!cdoc  drainage          GW-Abfluss
!cdoc  pEvSoil           max Evaporation bei Bodensaettigung
!cdoc  actEvap           aktuelle Bodenevaporation
!cdoc  rootdens          Wurzelverteilung 1/cm2 Bodenflaeche
!cdoc  van Genuchten
!cdoc  thetS             Porenvolumen m3/m3
!cdoc  thetR             min. Wassergehalt m3/m3
!cdoc  vgAlpha           1/Wassereintrittspunkt  cm-1
!cdoc  vgN               Porengroessenverteilungskoeffizient
!cdoc  kSat              gesaettigte Leitfaehigkeit cm/s
!cdoc  lambA-E           empirische Parameter fuer die BodenWaermeLeitfaehigkeit
!cdoc  1                 Grundwasser
!cdoc  SOM               SoilOrganicMatter g/g
!c---------------------------------- The End -----------------------------------
!cc==========================================================================
!c      function soilresp()
!cc
!cc     This function calculates the soil respiration by means of the
!cc     ARRHENIUS-Equation
!cc     sst: Temperature of Soil in Kelvin
!cc     r10: Scaling factor
!cc     ea: Activation Energy for Respiration
!cc--------------------------------------------------------------------------
!c      implicit none
!c
!c      integer nlay
!c      parameter (nlay=20)
!c
!c      common/const/pi,rgas,stef
!c      common/siteid/s_id
!c      common/vdriv/pyran,diff,dirn,tair,rh,sst,press,wind(nlay),ca
!c
!c      real rgas,stef
!c      real r10,ea,soilresp
!c      integer  s_id
!c      double precision pyran,diff,dirn,tair,rh,sst,press,wind,ca
!c
!c
!c      if (s_id.eq.10) then      !WEID
!cc---------------------------------------------------------
!cc         r10=2.957       ! residual, 2.column
!cc         ea=40000.
!cc---------------------------------------------------------
!c         r10=2.460        ! chamber measurements, 1.column
!c         ea=43430.
!c      elseif (s_id.eq.9) then    !THAR
!cc---------------------------------------------------------
!cc         r10=0.955       ! residual, 2.column
!cc         ea=55000.
!cc---------------------------------------------------------
!cc         r10=1.444        ! chamber measurements, 1.column
!cc         ea=43430.
!c         r10=1.570        ! adapted from tair
!c         ea=22360.
!c      elseif (s_id.eq.4) then    !HESS
!cc---------------------------------------------------------
!cc         r10=2.649       ! residual, 2.column
!cc         ea=46510.
!cc---------------------------------------------------------
!cc         r10=1.566        ! chamber measurements, 1.column
!cc         ea=73210.
!c         r10=2.539        ! adapted from tair
!c         ea=39260.
!c      elseif (s_id.eq.71) then   !VIBE
!cc---------------------------------------------------------
!cc         r10=2.273       ! residual, 2.column
!cc         ea=84590.
!cc---------------------------------------------------------
!cc         r10=2.568        ! chamber measurements, 1.column
!cc         ea=65510.
!c         r10=1.539        ! adapted from tair
!c         ea=97300.
!c      elseif (s_id.eq.8) then    !BRAS
!cc---------------------------------------------------------
!cc         r10=1.761       ! residual, 2.column
!cc         ea=21290.
!cc---------------------------------------------------------
!c         r10=1.081        ! chamber measurements, 1.column
!c         ea=73230.
!c      elseif (s_id.eq.6) then    !LOOB
!cc---------------------------------------------------------
!cc         r10=2.038       ! residual, 2.column
!cc         ea=58000.
!cc---------------------------------------------------------
!c         r10=3.400        ! chamber measurements, 1.column
!c         ea=43430.
!c      elseif (s_id.eq.5) then    !LILL
!cc---------------------------------------------------------
!cc         r10=3.059       ! residual, 2.column
!cc         ea=76720.
!cc---------------------------------------------------------
!c         r10=1.235        ! chamber measurements, 1.column
!c         ea=124970.
!c      elseif (s_id.eq.11) then   !ABER
!cc---------------------------------------------------------
!cc         r10=1.489       ! residual, 2.column
!cc         ea=180000.
!cc---------------------------------------------------------
!c         r10=1.081        ! chamber measurements, 1.column
!c         ea=73230.
!c      elseif (s_id.eq.14) then   !HYYT
!cc---------------------------------------------------------
!cc         r10=1.891       ! residual, 2.column
!cc         ea=48000.
!cc---------------------------------------------------------
!c         r10=1.000        ! chamber measurements, 1.column
!c         ea=43430.
!c      elseif (s_id.eq.72) then   !VIDO
!cc---------------------------------------------------------
!cc         r10=1.704       ! residual, 2.column
!cc         ea=60000.
!cc---------------------------------------------------------
!c         r10=1.284        ! chamber measurements, 1.column
!c         ea=65510.
!c      elseif (s_id.eq.50.or.s_id.eq.51) then   !GRAx
!cc         r10=2.46
!cc         ea=43430.
!c         r10=2.072        ! adapted from tair
!c         ea=35000.        ! Brache: 1.870 and 35000.
!cc---------------------------------------------------------
!c      elseif (s_id.eq.80) then   !MIXX, average from THAR and HESS
!cc---------------------------------------------------------
!c         r10=2.055        ! adapted from tair
!c         ea=30810.
!c      endif
!c
!cc     Unit umol m-2 s-1
!c      soilresp=r10*exp((ea*(sst-283.16))/(rgas*sst*283.16))
!c
!c      return
!c      end

!END MODULE PROXELMODULE
